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.2 ■ Abstract: We have generalized our "unified" model of evolutionary ecology by taking into 
• account the possible movements of the organisms from one "patch" to another within the 
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square lattice where each site corresponds to a distinct "patch" . A self-organizing hierarchical 
^ food web describes the prey-predator relations in the eco-system. The same species at different 
Q ' patches have identical food habits but differ from each other in their reproductive characteristic 
. features. By carrying out computer simulations up to 10^ time steps, we found that, depending 
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^ ' model turn out to be robust against migration of the organisms. 
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1 Introduction 



Direct experimental evidence in recent years [H |2l IHl IH El El El has established that significant 
evolutionary changes can occur in an ecosystem over ecologically relevant time scales. In other 
words, the dynamics of ecology and evolution are inseparable. Significant progress has been 
made over the last two years in developing detailed models that describe not only ecological 
phenomena on short periods of time but also evolutionary processes on longer time scales 
[3 El Uni HH Uni UH Uni Uni mi • some of these models describe the population dynamics of each 
species in terms of only a single dynamical variable, whereas some other more detailed models 
monitor the birth, ageing and death of individual organisms. A thorough comparison of these 
recent models will be given elsewhere |18j . 

In all the versions of our "unified" model jHl El QHl ^] so far we have assumed that the 
population of the prey as well as that of predators are uniformly distributed in space. The 
absence of the spatial degrees of freedom in such models of spatially-extended systems should 
be interpreted, using the language of statistical physics p[^, as a mean-field- like approximation. 
This situation is similar to a well-stirred chemical reaction where the spatial fiuctuations in the 
concentrations of the reactants and the products is negligibly small. On the other hand, spatial 
inhomogeneities in the eco-systems and migration of organisms from one eco-system to another 
are known to play crucial roles in evolutionary ecology OT] |22| • 

In recent years, the spatial inhomogeneity of the populations observed in real ecosystems 
have been captured by extending the Lotka-Volterra systems on discrete lattices where each of 
the lattice sites represents different spatial "patches" or "habitats" of the ecosystem [221 EH 
E3 ESI EH EH EH EOl EI] • However, almost all these investigations of "geographical" effects 
were restricted to only a single pair of species one of which is the predator while the other is 
the prey. 

For modelling the population dynamics of more than two species, one needs to know the food 
web which is a graphical way of describing the prey-predator relations, i.e., which species eats 
which one and which compete among themselves for the same food resources E3 IMl Ei] • 
More precisely, a food web is a directed graph where each node is labelled by a species' name 
and each directed link indicates the direction of flow of nutrient (i.e., from a prey to one of its 
predators) . 

A static (time-independent) food web may be a good approximation over a short period of 
time. But, a more realistic description, valid over longer periods of time, must take into account 
not only the adaptations of the species and their changing food habits, but also their extinction 
and creation of new species through speciation jSH]. These processes make the food web a 
slowly evolving graph. Such slow time evolution of the food webs were naturally incorporated 
in our earlier mean-fleld descriptions |Hl El UHl HI] • 

In this paper we extend the latest version of our "unifled" model ^1] of evolutionary ecology 
incorporating spatial inhomogeneities of the eco-system, i.e., spatial variations from one patch to 
another. Biologically motivated simulations on lattices, to take into account the geographical 
extent of an ecosystem, have a long tradition, e.g. for prey-predator relations ^37j, ageing 
[SH] . speciation [32]; their results sometimes [10] differ drastically compared with those of the 
homogeneous models that ignore the spatial extension. Such type of simulations are part of the 
widespread efforts to apply statistical physics methods to biology [H] or other flelds outside 
physics [121 • 
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Our model 13131101, in its latest and most realistic form [TT], allows each individual organism 
to give birth to M offspring at each time step, provided the age of the adult organism is 
above some minimum reproduction age Xrep (which varies from one species to another), with a 
probability that depends on the current age as well as X^ep- Similarly, each organism can die, 
for genetic reasons, with a probability (mortality rate) increasing exponentially with age, up to 
some species-dependent maximum age X^ax- Moreover, mutations happen at each time step, 
with the probability Pmuti which change X^ep and M randomly by ±1. 

In our model, we assumed a simple hierarchical food web where the species are arranged 
in trophic levels ^ (1 < ^ < (-max), with no more than m^~^ species in each level, i = 1 thus 
represents one species at the top of all food chains; typically, m = 2. Species in level i can 
prey on (or ignore) species on the immediately lower level i + 1. Moreover, the prey-predator 
relations are dynamic; because of random mutation, the food habits of the species change. 
Thus, an individual organism may fall victim to one of its predators. However, scarcity of food 
may also lead to starvation deaths among the predators. 

A species becomes extinct when all the corresponding population vanishes. We re-fill each 
ecological niche left empty by an extinct species, separately with probability Psp, by another 
species on the same (or a lower) level; independent random mutation of the parent and daughter 
species leads, eventually, to two different species. 

If the total biomass of the surviving individuals in the whole ecosystem is below the preas- 
signed maximum possible value, a new trophic level is created, with probability piev, i-e. imax 
is increased by one. The actual number of occupied levels can fluctuate with time depending 
on extinction and speciation. 

The maximum possible age of individual organisms, X^ax, for each species is assumed to 
depend only on the trophic level and fixed at 100/2*^^^^^/^. The minimum reproduction age Xrep 
is initially distributed randomly between one and this maximum age. We do not monitor the 
ageing and death of the the "bacteria" occupying the lowest level of the food web; instead, we 
assume a constant population of the "bacteria" throughout the evolution. 

The paper is organized as follows. In section 2 we study a simplified version of our "unified" 
model that provides a coarse-grained description of evolutionary ecology where the birth, ageing 
and death of the individual organisms do not appear explicitly. After studying the homogeneous 
version of this model we also investigate the effects of inhomogeneities arising from population 
migration. A few variants of this model are considered. In section 3 we extend the model 
formulating it in terms of the detailed "micro-dynamics" of individual organisms which can 
also migrate from one "patch" of the eco-system to another. The main results are summarized 
and conclusions are drawn in section 4. 

2 A coarse-grained description of evolutionary ecology 

In this section we do not monitor the "microscopic" dynamics of birth, ageing and death of the 
individuals explicitly. Instead, we describe the population dynamics of each species, in a coarse- 
grained manner, in terms of a single dynamical variable. In the first subsection below we ignore 
the spatial fluctuations while in the second subsection we take into account the fluctuations 
from one "patch" to another by modelling the eco-system as a lattice where each site represents 
a "patch". 
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2.1 Homogeneous approximation 
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Figure 1: Histogram for the lifetimes (part a) and minimum reproduction times (part b) of 
species, without lattice, with observation time t = 10^, lO'^, 10^, 10^, 10''' from right to left, 
summing over 10^, 10'^, 10^, 10, 1 samples, respectively. 
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Figure 2: Histogram for the lifetimes (part a) and minimum reproduction times (part b) of 
species, with diffusion on a 31 x 31 square lattice, with observation time t = 10^, 10^, 10^ from 
right to left, using one sample only. 

When we deal with species and no longer with individuals, we have to omit the individual 
births and genetic deaths and replace them by a fitness of the whole species; normalized to 
unity this fitness gives the survival probability of the species as a whole. Our fitness has two 
components, both taken from the birth and death probabilities of our previous model [TT] : 

a) High birth rates M decrease the life expectancy of the parent; thus a species is assumed 
to survive a time step with probability exp(— rM), where r is taken as 0.05 as in [II Ij . 

b) Having survived step a, a species is the fitter, the higher the number M of offspring per 
time step is and the lower the minimum age of reproduction Xrep is. Thus each species survives 
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Figure 3: Variation with time for the average of the maximal number of levels, with (x) and 
without (line) lattice (same simulations as in Figs.l to 4.). The + symbols show one run 
without lattice to much longer times. These data show the need to iterate up to 10^ time steps 
for equilibration. 



with probability [{Xmax — Xrep) M] / [{Xmax — Xmin)Mmax] whcrc Mmax IS the current maximum 
litter size M in the entire ecosystem, and Xmin the current minimum value of the Xrep in the 
whole ecosystem. In this way, the fittest species, and only this species, survives for sure, all 
others survive with the above probability only. 

[Another modification is required in the prey-predator deaths. The Verhulst deaths are now 
omitted, and instead of a reduction of the population by the amount given in eq(l) of JT], the 
population goes extinct with a probability given by this amount.] 

We took psp = 0.1, pn,ut = 0.001, piev = 0.0001, r = 0.05, C = 0.2 as in HH, and allowed 
the creation of a new level only if at most four species were alive above the lowest level. We 
averaged over t consecutive time steps, after making t/5 time steps to establish a reasonable 
equilibrium. 

We found that the number of food levels actually occupied at any one moment (between 1 
and 9) is normally distributed ("Gaussian" parabola on a semi - logarithmic plot; not shown). 
Thus for extremely long times t the maximum number of food levels should increase with the 
square- root of log(t). This extremely slow increase is compatible with one simulation (Figs. 3 
and 7) for t = 10^. Biologically, longer times do not make much sense; 100 million years ago 
we did not have the genus homo on earth. 

The histograms for the lifetimes of species, and for the minimum reproduction times Xrep, 
Fig.l, both show that t = 10^ and 10^ are too short, t = 10^ is intermediate between short- 
and long-time behaviour, and t = 10^ and 10^ give barely distinguishable long-time behaviour. 
For this comparison we summed up 10^/t samples so that the total number of measurements 
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is the same in all five cases. 

Fig. la also shows an exponential lifetime histogram for long t, as opposed to the more 
complicated histogram in jHUHlElE] and the simple power law often claimed for reality 
Thus the treatment of individuals, not the other complications, seems to be responsible for the 
complex lifetime distributions in these earlier models. Fig. lb, on the other hand, is not much 
different from Fig.3 in [llj. The average litter size M fluctuates strongly, increasing overall 
perhaps logarithmically with the number of time steps (not shown). 

2.2 Spatial inhomogeneity 

With a whole ecosystem of several food layers on every site of a lattice, computational re- 
quirements increase drastically. To reduce memory requirements we thus assume that the 
prey-predator relations are the same on each lattice site; only the other species-dependent 
properties like Xrep and M, and the number i^ax of food levels, can change from site to site. 
Mutations thus act on the latter quantities differently for different sites, while mutations on 
prey-predator relations act simultaneously on the whole lattice, for the given niche (= possible 
species in ecosystem). If the expansion of one species to another site, where it later undergoes 
different mutations, is called speciation (i.e. parapatric speciation), then the previously dis- 
cussed speciation events which also change the prey-predator relations, are changes of a higher 
taxonomical level, like changes of genus. Alternatively, if we call the changes of the preceding 
section a speciation, the changes due to expansion onto two lattice sites allows the formation 
of sub-species ( "races" ) of the same species, sharing the same predator-prey relations. 

If we only allow for different ecosystems on different lattices sites, the results are roughly 
equivalent to averaging over many separate runs of a single site, as we did in the previous 
section. Only the fluctuations are reduced due to the agreement for prey-predator relations. 
More interesting and realistic is motion of species across the lattice. Thus at every time step 
each non-extinct species with probability 1/2 stays where it is, and otherwise selects randomly 
one of its four neighbour sites on the square lattice (using periodic boundary conditions). If 
that site is occupied at the corresponding niche (i.e. at the same level and on the same node) 
then again nothing happens; if the selected neighbour site is not already occupied for that 
niche, then the old species occupies this place but also stays at its old site. Thus we have 
an expansion, not a mere displacement, of the species. From then on, both sites undergo the 
usual mutations and thus slowly the two populations drift apart genetically to become different 
(sub-)species. 

Apart from somewhat shorter lifetimes, our results in Fig. 2 look similar to Figs. 1: Simple 
exponential decay for the distributions of lifetimes and minimum reproduction ages, for long 
enough times. More interesting is Fig.3 with the maximal number of levels, imax, averaged over 
all lattice sites. As a function of time, this quantity at first increases roughly linearly in time 
from its initial value 3. After about 10^ time steps, a much slower increase is observed, leading 
to a plateau near 8. Diffusion barely changes the lattice results without diffusion, and even 
without a lattice the results are about the same except that the single run without a lattice 
shows stepwise increase of imaxit)- 



6 




Figure 4: As Fig.l, but with fitness tested every tenth time step, and r ten times smaller. The 
two longest times are shown by small dots in the left part. The straight line has the slope -2 
sometimes claimed for real extinctions. 




Figure 5: With diffusion on a 31 x 31 square lattice, / = 0.01, r = 0.5. 
2.3 Variants 

In the above simulations, the lowest level is by construction always fully occupied, but only a 
few species exist in the higher levels. We avoided the extinction of so many species by reducing 
the "birth punishment" factor r from 0.05 to 0.005, and by making the fitness test only with 
a probability of /, that means, on average only at every (l//)-th time step the less fit species 
may go extinct. Now figure 4 for / = 0.1 shows for the lifetime distribution a more complicated 
behaviour, similar to [HliniCniEl; without a lattice. For diffusion on a lattice, the lifetimes 
may become log-normally distributed for these parameters. Fig. 5 shows a more interesting 
distribution for / = 0.01, r = 0.5. Simulations with different system size give similar results. 
Fig. 6. In real palaeontology, of course, it would be difficult to measure lifetime distributions 
over ten decades with a single criterion, as is done here and in Fig. 2. Fig. 7 suggests that the 
maximum number of layers increases very weakly; again times varying over nine decades are 
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L = 1, 2, 5, 31, 61 (bottom to top); t=10"5 



L= 1, 2, 5, 31, 61 (bottom to top); t=10"5 
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Figure 6: Parameters as in Fig. 2 (Fig. 6a, r = 0.05, without fitness test), and 4 (Fig.Gb, 
r = 0.005,/ = 0.1) for various lattice sizes L at t = 10^. Since only one run was made, the 
resulting histograms increase with increasing L. 



difficult to observe palaeontologially. 



3 Inhomogeneous eco-system: a individual-based descrip- 
tion 

In this section we report results of the most detailed description of the eco-system. We not 
only incorporate explicitly the birth, ageing and death of individual organisms, but also take 
into account the spatial inhomogeneities of the populations in terms of a lattice model that is 
very similar to the one considered in the preceding section. 

Thus we assumed that each individual either ages by one time unit for each time step, or 
it dies. Besides the deaths as prey, unchanged from the previous section, the above selection 
according to fitness and litter size is replaced by a Gompertz-type death probability as a function 
of the age a of the individual, and depending also on the maximum age X^ax and litter size M 
of the whole species: 

Pdeath = exp[(a - X.raax)r/M] 

where r is the same free parameter, e.g. 0.05, as in the earlier sections. (For ages below the 
minimum reproduction age Xrep the death probability is assumed to be age-independent, with 
a replaced by Xrep in the above equation.) 

The survivors give birth to M offspring at each time step with probability 

Pbirth (^maa; '^)/ i^Xmax X^ep) 

if their age a is above the minimum reproduction age Xrep- The just-born babies die with a 
Verhulst probability jlHI ^—n/umax where n is the actual current population of the given species 
on the given lattice site, and Umax, taken here as 100, is its maximum possible population. 

A population may expand into a neighbouring lattice site, if the population there is zero 
for the same species on the same level. The population stays the same at the old site, whereas 
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Figure 7: Variation of imax with time up to nearly 10^ time steps for the parameters of Fig.4, 
showing the less than logarithmic increase at long times (big symbols). The small dots below 
the big symbols are the actual numbers of occupied levels in the same simulations. The higher 
dashed line gives imax for ten times higher probabilities of mutation and level formation. 

at the new site, as in our standard sympatric speciation assumption, the population is taken 
randomly between one and the population at the old site; also the age of the emigrants starts 
at zero. 

Fig. 8 shows lifetimes and minimum age of reproduction with the individuals on a 31 x 31 
lattice; smaller lattices gave similar curves. To speed up equilibration, pmut = 0.01, piev = 0.001 
were taken ten times higher than before. To prevent the number of levels to increase towards 
10, a new level could be created only if the total population on that lattice site was at most 
'^TT'max = 200. We See that roughly but not precisely the lifetimes follow an inverse-squared 
power law (straight line in Fig.Sa). 

4 Discussion 

We have extended our earher "unified" models of evolutionary ecology jHl El UDl HI] by taking 
into account the spatial inhomogeneities arising from migration of populations from one patch 
to another. Besides, individuals and species, we also allowed evolution of genera (i.e., sets of 
species with the same food habits). By carrying out extensive computer simulations, we find 
that, depending on the details of the parameter set, the distribution of the lifetimes of the 
species may be simple exponential or, as before, a combination of power laws. Other aspects 
like, for example, the self-organization of the minimum reproduction ages, are less affected by 
this generalization. 
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Figure 8: Distribution of lifetimes and minimum reproduction ages with individuals on a 13 
lattice at t = 10''^, iC, 10^ 10^. The thin curves (or dots) refer tot = 10^ for 5 x 5 (lower data) 
and 31 X 31 (higher data) for the lifetimes. 
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